/******************************************************************************
 * bn.c - big number math implementation
 *
 * Copyright (c) 2004 by Juergen Buchmueller <pullmoll@stop1984.com>
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software Foundation,
 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA
 *
 *	$Id: bignum.c,v 1.19 2006/02/21 16:14:26 pullmoll Exp $
 ******************************************************************************/
#include <math.h>
#include "bignum.h"

#if	BN_TEST
#define	DBG(x)	info x
#define	bn_snprintf snprintf
#else
#define	DBG(x)
#define	bn_snprintf pm_snprintf
#endif

#ifndef	ASM_X86
#ifdef	X86
#define	ASM_X86	1
#endif
#endif

/**
 * @brief The upper two bits of a limb (used in RSA functions)
 */
#define	DIGMSB2 (DIGMSB|(DIGMSB>>1))

/**
 * @brief A table of the first 1999 odd primes
 */
static limb_t sm_primes[] = {
	   3,   5,   7,  11,  13,  17,  19,  23,  29,  31,
	  37,  41,  43,  47,  53,  59,  61,  67,  71,  73,
	  79,  83,  89,  97, 101, 103, 107, 109, 113, 127,
	 131, 137, 139, 149, 151, 157, 163, 167, 173, 179,
	 181, 191, 193, 197, 199, 211, 223, 227, 229, 233,
	 239, 241, 251, 257, 263, 269, 271, 277, 281, 283,
	 293, 307, 311, 313, 317, 331, 337, 347, 349, 353,
	 359, 367, 373, 379, 383, 389, 397, 401, 409, 419,
	 421, 431, 433, 439, 443, 449, 457, 461, 463, 467,
	 479, 487, 491, 499, 503, 509, 521, 523, 541, 547,
	 557, 563, 569, 571, 577, 587, 593, 599, 601, 607,
	 613, 617, 619, 631, 641, 643, 647, 653, 659, 661,
	 673, 677, 683, 691, 701, 709, 719, 727, 733, 739,
	 743, 751, 757, 761, 769, 773, 787, 797, 809, 811,
	 821, 823, 827, 829, 839, 853, 857, 859, 863, 877,
	 881, 883, 887, 907, 911, 919, 929, 937, 941, 947,
	 953, 967, 971, 977, 983, 991, 997,1009,1013,1019,
	1021,1031,1033,1039,1049,1051,1061,1063,1069,1087,
	1091,1093,1097,1103,1109,1117,1123,1129,1151,1153,
	1163,1171,1181,1187,1193,1201,1213,1217,1223,1229,
	1231,1237,1249,1259,1277,1279,1283,1289,1291,1297,
	1301,1303,1307,1319,1321,1327,1361,1367,1373,1381,
	1399,1409,1423,1427,1429,1433,1439,1447,1451,1453,
	1459,1471,1481,1483,1487,1489,1493,1499,1511,1523,
	1531,1543,1549,1553,1559,1567,1571,1579,1583,1597,
	1601,1607,1609,1613,1619,1621,1627,1637,1657,1663,
	1667,1669,1693,1697,1699,1709,1721,1723,1733,1741,
	1747,1753,1759,1777,1783,1787,1789,1801,1811,1823,
	1831,1847,1861,1867,1871,1873,1877,1879,1889,1901,
	1907,1913,1931,1933,1949,1951,1973,1979,1987,1993,
	1997,1999,2003,2011,2017,2027,2029,2039,2053,2063,
	2069,2081,2083,2087,2089,2099,2111,2113,2129,2131,
	2137,2141,2143,2153,2161,2179,2203,2207,2213,2221,
	2237,2239,2243,2251,2267,2269,2273,2281,2287,2293,
	2297,2309,2311,2333,2339,2341,2347,2351,2357,2371,
	2377,2381,2383,2389,2393,2399,2411,2417,2423,2437,
	2441,2447,2459,2467,2473,2477,2503,2521,2531,2539,
	2543,2549,2551,2557,2579,2591,2593,2609,2617,2621,
	2633,2647,2657,2659,2663,2671,2677,2683,2687,2689,
	2693,2699,2707,2711,2713,2719,2729,2731,2741,2749,
	2753,2767,2777,2789,2791,2797,2801,2803,2819,2833,
	2837,2843,2851,2857,2861,2879,2887,2897,2903,2909,
	2917,2927,2939,2953,2957,2963,2969,2971,2999,3001,
	3011,3019,3023,3037,3041,3049,3061,3067,3079,3083,
	3089,3109,3119,3121,3137,3163,3167,3169,3181,3187,
	3191,3203,3209,3217,3221,3229,3251,3253,3257,3259,
	3271,3299,3301,3307,3313,3319,3323,3329,3331,3343,
	3347,3359,3361,3371,3373,3389,3391,3407,3413,3433,
	3449,3457,3461,3463,3467,3469,3491,3499,3511,3517,
	3527,3529,3533,3539,3541,3547,3557,3559,3571,3581,
	3583,3593,3607,3613,3617,3623,3631,3637,3643,3659,
	3671,3673,3677,3691,3697,3701,3709,3719,3727,3733,
	3739,3761,3767,3769,3779,3793,3797,3803,3821,3823,
	3833,3847,3851,3853,3863,3877,3881,3889,3907,3911,
	3917,3919,3923,3929,3931,3943,3947,3967,3989,4001,
	4003,4007,4013,4019,4021,4027,4049,4051,4057,4073,
	4079,4091,4093,4099,4111,4127,4129,4133,4139,4153,
	4157,4159,4177,4201,4211,4217,4219,4229,4231,4241,
	4243,4253,4259,4261,4271,4273,4283,4289,4297,4327,
	4337,4339,4349,4357,4363,4373,4391,4397,4409,4421,
	4423,4441,4447,4451,4457,4463,4481,4483,4493,4507,
	4513,4517,4519,4523,4547,4549,4561,4567,4583,4591,
	4597,4603,4621,4637,4639,4643,4649,4651,4657,4663,
	4673,4679,4691,4703,4721,4723,4729,4733,4751,4759,
	4783,4787,4789,4793,4799,4801,4813,4817,4831,4861,
	4871,4877,4889,4903,4909,4919,4931,4933,4937,4943,
	4951,4957,4967,4969,4973,4987,4993,4999,5003,5009,
	5011,5021,5023,5039,5051,5059,5077,5081,5087,5099,
	5101,5107,5113,5119,5147,5153,5167,5171,5179,5189,
	5197,5209,5227,5231,5233,5237,5261,5273,5279,5281,
	5297,5303,5309,5323,5333,5347,5351,5381,5387,5393,
	5399,5407,5413,5417,5419,5431,5437,5441,5443,5449,
	5471,5477,5479,5483,5501,5503,5507,5519,5521,5527,
	5531,5557,5563,5569,5573,5581,5591,5623,5639,5641,
	5647,5651,5653,5657,5659,5669,5683,5689,5693,5701,
	5711,5717,5737,5741,5743,5749,5779,5783,5791,5801,
	5807,5813,5821,5827,5839,5843,5849,5851,5857,5861,
	5867,5869,5879,5881,5897,5903,5923,5927,5939,5953,
	5981,5987,6007,6011,6029,6037,6043,6047,6053,6067,
	6073,6079,6089,6091,6101,6113,6121,6131,6133,6143,
	6151,6163,6173,6197,6199,6203,6211,6217,6221,6229,
	6247,6257,6263,6269,6271,6277,6287,6299,6301,6311,
	6317,6323,6329,6337,6343,6353,6359,6361,6367,6373,
	6379,6389,6397,6421,6427,6449,6451,6469,6473,6481,
	6491,6521,6529,6547,6551,6553,6563,6569,6571,6577,
	6581,6599,6607,6619,6637,6653,6659,6661,6673,6679,
	6689,6691,6701,6703,6709,6719,6733,6737,6761,6763,
	6779,6781,6791,6793,6803,6823,6827,6829,6833,6841,
	6857,6863,6869,6871,6883,6899,6907,6911,6917,6947,
	6949,6959,6961,6967,6971,6977,6983,6991,6997,7001,
	7013,7019,7027,7039,7043,7057,7069,7079,7103,7109,
	7121,7127,7129,7151,7159,7177,7187,7193,7207,7211,
	7213,7219,7229,7237,7243,7247,7253,7283,7297,7307,
	7309,7321,7331,7333,7349,7351,7369,7393,7411,7417,
	7433,7451,7457,7459,7477,7481,7487,7489,7499,7507,
	7517,7523,7529,7537,7541,7547,7549,7559,7561,7573,
	7577,7583,7589,7591,7603,7607,7621,7639,7643,7649,
	7669,7673,7681,7687,7691,7699,7703,7717,7723,7727,
	7741,7753,7757,7759,7789,7793,7817,7823,7829,7841,
	7853,7867,7873,7877,7879,7883,7901,7907,7919
};

/**
 * @brief number of small primes to test as divisors
 */
#define	SM_PRIMES	(sizeof(sm_primes)/sizeof(sm_primes[0]))

/**
 * @brief number of small primes to test before Rabin-Miller
 */
#define	SM_FIRST	512

/**
 * @brief Return hex representation of a big number
 *
 * Returns the hex representation of a[],
 * where a is a big number integer with nlimb limbs.
 *
 * @param a pointer to an array of limbs
 * @param nlimb number of limbs in the array
 *
 * @result string containing the hex representation of a
 */
const char *bn2x(limb_t *a, cnt_t nlimb)
{
	static cnt_t which = 0;
	static char *xbuff[8] = {
		NULL, NULL, NULL, NULL,
		NULL, NULL, NULL, NULL
	};
	char *dst;
	cnt_t size;
	cnt_t n = nlimb;

	if (0 == n)
		return "0";

	which = (which + 1) % 8;
	size = 8 * n + 1;
	if (NULL != xbuff[which])
		free(xbuff[which]);
	dst = xbuff[which] = calloc(size, sizeof(char));
	if (NULL == dst)
		return "memory error";
	while (n-- > 0) {
		dst += bn_snprintf(dst, size, "%08x", a[n]);
		size -= 8;
	}
	return xbuff[which];
}

/**
 * @brief Return decimal representation of a big number
 *
 * Returns the decimal representation of a[],
 * where a is a big number integer with nlimb limbs.
 *
 * @param a pointer to an array of limbs
 * @param nlimb number of limbs in the array
 *
 * @result string containing the decimal representation of a
 */
const char *bn2d(limb_t *a, cnt_t nlimb)
{
	static cnt_t which = 0;
	static char *dbuff[8] = {
		NULL, NULL, NULL, NULL,
		NULL, NULL, NULL, NULL
	};
	static limb_t v[BN_MAXSIZE];
	limb_t r;
	char *dst;
	cnt_t size;
	cnt_t n = bn_sizeof(a, nlimb);

	if (0 == n)
		return "0";

	bn_copy(v, a, n);
	which = (which + 1) % 8;
	size = 12 * n + 1;
	if (NULL != dbuff[which])
		free(dbuff[which]);
	dst = dbuff[which] = calloc(size, sizeof(char));
	if (NULL == dst)
		return "memory error";
	size--;
	while (0 != bn_cmp_limb(v, 0, n)) {
		r = bn_div_limb(v, v, 10, n);
		dst[--size] = '0' + r;
	}
	return &dst[size];
}

/**
 * @brief Return decimal representation of a big number pair
 *
 * Returns the decimal representation of a[].b[],
 * where a is a big number integer with alimb limbs,
 * and b is a multiprecision fixed fraction with blimb limbs.
 *
 * @param a pointer to an array of limbs
 * @param alimb number of limbs in the a array
 * @param b pointer to an array of limbs
 * @param blimb number of limbs in the b array
 *
 * @result string containing the decimal representation of a.b
 */
const char *bn2f(limb_t *a, cnt_t alimb, limb_t *b, cnt_t blimb)
{
	static cnt_t which = 0;
	static char *dbuff[8] = {
		NULL, NULL, NULL, NULL,
		NULL, NULL, NULL, NULL
	};
	static limb_t v[BN_MAXSIZE];
	static limb_t w[BN_MAXSIZE];
	limb_t r;
	char *dst;
	cnt_t size;

	bn_copy(v, a, alimb);
	bn_copy(w, b, blimb);

	which = (which + 1) % 8;
	size = 12 * (alimb + blimb) + 1 + 1;
	if (NULL != dbuff[which])
		free(dbuff[which]);
	dst = dbuff[which] = calloc(size, sizeof(char));
	if (NULL == dst)
		return "memory error";
	size = 12 * alimb;
	while (0 != bn_cmp_limb(w, 0, blimb) && size < 12 * (alimb + blimb)) {
		r = bn_mul_limb(w, w, 10, blimb);
		dst[size++] = '0' + r;
	}

	size = 12 * alimb;
	dst[size] = '.';
	while (0 != bn_cmp_limb(v, 0, alimb) && size > 0) {
		r = bn_div_limb(v, v, 10, alimb);
		dst[--size] = '0' + r;
	}

	return &dst[size];
}

/**
 * @brief Return binary representation of a big number
 *
 * Returns the binary representation of a[],
 * where a is a big number integer with nlimb limbs.
 *
 * @param a pointer to an array of limbs
 * @param nlimb number of limbs in the array
 *
 * @result string containing the binary representation of a
 */
const char *bn2b(limb_t *a, cnt_t nlimb)
{
	static cnt_t which = 0;
	static char *bbuff[8] = {
		NULL, NULL, NULL, NULL,
		NULL, NULL, NULL, NULL
	};
	limb_t r;
	char *dst;
	cnt_t size;
	cnt_t n = bn_sizeof(a, nlimb);

	if (0 == n)
		return "0";

	which = (which + 1) % 8;
	size = LIMBBITS * n + 1;
	if (NULL != bbuff[which])
		free(bbuff[which]);
	dst = bbuff[which] = calloc(size, sizeof(char));
	if (NULL == dst)
		return "memory error";
	n = 0;
	size--;
	while (size-- > 0) {
		r = (a[n/LIMBBITS] >> (n%LIMBBITS)) & 1;
		n++;
		dst[size] = '0' + r;
	}
	return &dst[size];
}

/**
 * @brief Zero an array of limbs
 *
 * Sets a[] = 0
 * where a is a big number integer of nlimb limbs.
 *
 * @param a pointer to an array of limbs
 * @param nlimb number of limbs in the array
 *
 */
void bn_zero(limb_t a[], cnt_t nlimb)
{
	memset(a, 0, nlimb * sizeof(limb_t));
}

/**
 * @brief Set an array of limbs to a single limb value
 *
 * Sets a[] = d
 * where a is a big number integer of nlimb limbs,
 * and d is a single limb
 *
 * @param a pointer to an array of limbs to set
 * @param d limb value to set a to
 * @param nlimb number of limbs in the array
 *
 */
void bn_limb(limb_t a[], limb_t d, cnt_t nlimb)
{
	memset(a, 0, nlimb * sizeof(limb_t));
	a[0] = d;
}

/**
 * @brief Copy an array of limbs
 *
 * Sets a[] = b[]
 * where a and b are a big number integers of nlimb limbs
 *
 * @param a pointer to an array of limbs (destination)
 * @param b pointer to an array of limbs (source)
 * @param nlimb number of limbs in the arrays
 */
void bn_copy(limb_t a[], limb_t b[], cnt_t nlimb)
{
	memcpy(a, b, nlimb * sizeof(limb_t));
}

/**
 * @brief Return significant size of a big number
 *
 * Returns size of significant limbs in a[]
 * i.e. searches for the first non-zero limb from
 * nlimb-1 downto 0.
 *
 * @param a pointer to an array of limbs (candidate)
 * @param nlimb number of limbs in the arrays
 *
 * @result number of significant limbs in a
 */
cnt_t bn_sizeof(limb_t a[], cnt_t nlimb)
{
	while (nlimb-- > 0)
		if (0 != a[nlimb])
			return ++nlimb;
	return 0;
}


/**
 * @brief Return sign of a bignum minus a limb
 *
 * Returns the sign of (a[] - b)
 * where a is a big number integer of nlimb limbs,
 * and b is a single limb
 +
 * @param a pointer to an array of limbs (minuend)
 * @param b a single limb (subtrahend)
 * @param nlimb number of limbs in the array a
 *
 * @result sign of the comparison: -1 a<b, 0 a=b, +1 a>b
 */
int bn_cmp_limb(limb_t a[], limb_t b, cnt_t nlimb)
{
	if (0 == nlimb)
		return 0;

	while (nlimb-- > 1)
		if (0 != a[nlimb])
			return +1;
	if (a[0] < b)
		return -1;
	if (a[0] > b)
		return +1;
	return 0;
}

/**
 * @brief Return sign of bignum a minus bignum b
 *
 * Returns the sign of (a[] - b[])
 * where a and b are a big number integers of nlimb limbs
 *
 * @param a pointer to an array of limbs (minuend)
 * @param b pointer to an array of limbs (subtrahend)
 * @param nlimb number of limbs in the arrays
 *
 * @result sign of the comparison: -1 a<b, 0 a=b, +1 a>b
 */
int bn_cmp(limb_t a[], limb_t b[], cnt_t nlimb)
{
	if (0 == nlimb)
		return 0;

	while (nlimb-- > 0) {
		if (a[nlimb] > b[nlimb])
			return +1;	/* GT */
		if (a[nlimb] < b[nlimb])
			return -1;	/* LT */
	}

	return 0;	/* EQ */
}

/**
 * @brief Single limb is even test
 *
 * Returns 1 if a is even, else 0
 * where a is a single limb
 *
 * @param a a single limb
 *
 * @result zero if a is odd, 1 if a is even
 */
int sl_iseven(limb_t a)
{
	return (a & 1) ? 0 : 1;
}

/**
 * @brief bignum is even test
 *
 * Returns 1 if a[] is even, else 0
 * where a is a big number integer of nlimb limbs
 * Note: a zero limb big number integer is even!
 *
 * @param a pointer to an array of limbs
 * @param nlimb number of limbs in the arrays
 *
 * @result zero if a is odd, 1 if a is even
 */
int bn_iseven(limb_t *a, cnt_t nlimb)
{
	if (0 == nlimb)
		return 1;
	return (a[0] & 1) ? 0 : 1;
}

/**
 * @brief Add a single limb to a bignum
 *
 * Computes w[] = u[] + v
 * where w, u are big number integers of nlimb lims each,
 * and v is a single limb.
 * Returns carry if the addition overflows.
 *
 * Ref: Derived from Knuth Algorithm A.
 *
 * @param w pointer to an array of limbs receiving result
 * @param u pointer to an array of limbs (addend 1)
 * @param v a single limb
 * @param nlimb number of limbs in the arrays w and u
 *
 * @result The carry status of the addition
 */
limb_t bn_add_limb(limb_t w[], limb_t u[], limb_t v, cnt_t nlimb)
{
	limb_t carry;
	cnt_t j;

	/* Copy u to w, so we can bail out if no borrow is left */
	if (w != u)
		bn_copy(w, u, nlimb);

	/* Add v to first limb of u */
	w[0] += v;
	carry = (w[0] < v ? 1 : 0);

	/* Add carry to subsequent limbs */
	for (j = 1; 0 != carry && j < nlimb; j++) {
		w[j] += carry;
		carry = (w[j] < carry ? 1 : 0);
	}
	return carry;
}


/**
 * @brief Subtract a single limb from a bignum
 *
 * Computes w[] = u[] - v
 * where w, u are big number integers of nlimb limbs each,
 * and v is a single limb.
 * Returns borrow (0 if u >= v, or 1 if v > u).
 * 
 * Ref: Derived from Knuth Algorithm S.
 *
 * @param w pointer to an array of limbs receiving the result
 * @param u pointer to an array of limbs (minuend)
 * @param v single limb (subtrahend)
 * @param nlimb number of limbs in the arrays
 *
 * @result borrow of the subtraction (0 if u >= v, 1 if u < v)
 */
limb_t bn_sub_limb(limb_t w[], limb_t u[], limb_t v, cnt_t nlimb)
{
	limb_t borrow;
	cnt_t j;

	/* Copy u to w, so we can bail out if no borrow is left */
	if (w != u)
		bn_copy(w, u, nlimb);

	/* Subtract v from first limb of u */
	w[0] -= v;
	borrow = (w[0] > ~v ? 1 : 0);

	/* Subtract borrow from subsequent limbs */
	for (j = 1; 0 != borrow && j < nlimb; j++) {
		w[j] -= borrow;
		borrow = (w[j] > ~borrow ? 1 : 0);
	}

	return borrow;
}

/**
 * @brief Divide a bignum by a single limb
 *
 * Computes quotient q[] = u[] / v
 * and returns remainder r = u[] % v
 * where q, u are big number integers of nlimb limbs each,
 * and v is a single limb.
 * 
 * Makes no assumptions about normalisation.
 * 
 * Ref: Knuth Vol 2 Ch 4.3.1 Exercise 16 p625
 *
 * @param q pointer to an array of limbs receiving the quotient
 * @param u pointer to an array of limbs (dividend)
 * @param v single limb (divisor)
 * @param nlimb number of limbs in the arrays
 *
 * @result single limb remainder of the division (modulo)
 */
limb_t bn_div_limb(limb_t q[], limb_t u[], limb_t v, cnt_t nlimb)
{
	cnt_t j;
	limb_t t[2], r;
	cnt_t shift;

	if (0 == nlimb)
		return 0;
	if (0 == v)
		return LIMBMASK;	/* Divide by zero error */

	/*
	 * Normalize first:
	 * qequires high bit of V to be set,
	 * so find most significant by shifting
	 * until DIGMSB is set.
	 */
	for (shift = 0; 0 == (v & DIGMSB); shift++)
		v <<= 1;
	r = bn_shl(q, u, shift, nlimb);
	
	j = nlimb;
	while (j-- > 0) {
		t[0] = q[j];
		t[1] = r;
		sl_div(&q[j], &r, t, v);
	}

	/* Unnormalize */
	r >>= shift;
	return r;
}

/**
 * @brief Modulo a bignum by a single limb
 *
 * Computes remainder (modulo) r = u[] mod v
 * Computes r = u[] mod v
 * where u is a big number integer of nlimb
 * and r, v are single precision limbs
 * 
 * Use remainder from divide function.
 *
 * @param u pointer to an array of limbs (dividend)
 * @param v single limb (divisor)
 * @param nlimb number of limbs in the arrays
 *
 * @result single limb remainder of the division (modulo)
 */
limb_t bn_mod_limb(limb_t u[], limb_t v, cnt_t nlimb)
{
	static limb_t q[2*BN_MAXSIZE];
	limb_t r;

	r = bn_div_limb(q, u, v, nlimb);

	bn_zero(q, nlimb);
	return r;
}

/**
 * @brief Multiply a bignum by a single limb
 *
 * Computes product w[] = u[] * v
 * Returns overflow k
 * where w, u are big number integers of nlimb each
 * and v is a single limb
 *
 * @param w pointer to an array of limbs to receive the result
 * @param u pointer to an array of limbs (factor)
 * @param v single limb (other factor)
 * @param nlimb number of limbs in the arrays
 *
 * @result zero if no overflow, else overflow (value of w[nlimb])
 */
limb_t bn_mul_limb(limb_t w[], limb_t u[], limb_t v, cnt_t nlimb)
{
	limb_t t[2];
	limb_t carry;
	cnt_t j;

	if (0 == v) {
		bn_zero(w, nlimb);
		return 0;
	}

	for (j = 0, carry = 0; j < nlimb; j++) {
		sl_mul(t, u[j], v);
		w[j] = t[0] + carry;
		carry = t[1] + (w[j] < carry ? 1 : 0);
	}

	return carry;
}

#if	HAVE_U64
/**
 * @brief Computes quotient and remainder of 64 bit / 32 bit
 *
 * Computes quotient q = u[] / v, remainder r = u[] mod v
 * where u[] is a double limb.
 *
 * With native support for double limb division
 *
 * @param q pointer to the limb to receive the quotient
 * @param r pointer to the limb to receive the remainder
 * @param u pointer to an array of two limbs
 * @param v single limb divisor
 *
 * @result zero on success
 */
limb_t sl_div(limb_t *q, limb_t *r, limb_t u[2], limb_t v)
{
#if	ASM_X86
	limb_t qq;
	limb_t rr;

	if (0 == v)
		/* division by zero */
		return LIMBMASK;
	asm volatile(
		"divl	%4"
		: "=a"(qq), "=d"(rr)
		: "a"(u[0]), "d"(u[1]), "g"(v));
	*q = qq;
	*r = rr;
#else
	dlimb_t dd;

	if (0 == v)
		/* division by zero */
		return LIMBMASK;
	dd = ((dlimb_t)u[1] << LIMBBITS) | u[0];
	*q = dd / v;
	*r = dd % v;
#endif
	return 0;
}

#else

#define B (HALFMASK + 1)

/**
 * @brief Computes quotient and remainder of 64 bit / 32 bit
 *
 * Computes quotient q = u / v, remainder r = u mod v
 * where u is a double limb
 * and q, v, r are single precision limbs.
 * Returns high limb of quotient (max value is 1)
 * Assumes normalized such that v1 >= b/2
 * where b is size of HALF_DIGIT
 * i.e. the most significant bit of v should be one
 * 
 * In terms of half-limbs in Knuth notation:
 *   (q2q1q0) = (u4u3u2u1u0) / (v1v0)
 *   (r1r0) = (u4u3u2u1u0) % (v1v0)
 * for m = 2, n = 2 where u4 = 0
 *
 * We set q = (q1q0) and return q2 as "overflow'
 * Returned q2 is either 0 or 1.
 *
 * @param q pointer to the limb to receive the quotient
 * @param r pointer to the limb to receive the remainder
 * @param u pointer to an array of two limbs
 * @param v single limb divisor
 *
 * @result zero on success
 */
limb_t sl_div(limb_t *q, limb_t *r, limb_t u[2], limb_t v)
{
	limb_t quot;
	limb_t rem;
	limb_t ul;
	limb_t uh;
	limb_t p0;
	limb_t p1;
	limb_t v0;
	limb_t v1;
	limb_t u0;
	limb_t u1;
	limb_t u2;
	limb_t u3;
	limb_t borrow;
	limb_t q1;
	limb_t q2;
	limb_t s;
	limb_t t;

	/* Check for normalisation */
	if (0 == (v & DIGMSB)) {
		*q = *r = 0;
		return LIMBMASK;
	}
	
	/* Split up into half-limbs */
	v0 = LSH(v);
	v1 = MSH(v);
	u0 = LSH(u[0]);
	u1 = MSH(u[0]);
	u2 = LSH(u[1]);
	u3 = MSH(u[1]);

	/* Do three rounds of Knuth Algorithm D Vol 2 p272 */

	/*
	 * ROUND 1 calculate q2:
	 * estimate quot = (u4u3)/v1 = 0 or 1,
	 * then set (u4u3u2) -= quot*(v1v0) where u4 = 0.
	 */
	quot = u3 / v1;
	if (quot > 0) {
		rem = u3 - quot * v1;
		t = SHL(rem) | u2;
		if (quot * v0 > t)
			quot--;
	}
	uh = 0;		/* (u4) */
	ul = u[1];	/* (u3u2) */
	if (quot > 0) {
		/* (u4u3u2) -= quot*(v1v0) where u4 = 0 */
		p0 = quot * v0;
		p1 = quot * v1;
		s = p0 + SHL(p1);
		ul -= s;
		borrow = (ul > ~s ? 1 : 0);
		uh -= MSH(p1) - borrow;

		if (0 != MSH(uh)) {
			/* add back */
			quot--;
			ul += v;
			uh = 0;
		}
	}
	q2 = quot;

	/*
	 * ROUND 2 calculate q1:
	 * estimate quot = (u3u2) / v1,
	 * then set (u3u2u1) -= quot*(v1v0)
	 */
	t = ul;
	quot = t / v1;
	rem = t - quot * v1;
	/* Test on v0 */
	t = SHL(rem) | u1;
	if (B == quot || (quot * v0) > t) {
		quot--;
		rem += v1;
		t = SHL(rem) | u1;
		if (rem < B && (quot * v0) > t)
			quot--;
	}

	/*
	 * multiply and subtract:
	 * (u3u2u1)' = (u3u2u1) - quot*(v1v0)	
	 */
	uh = MSH(ul);	/* (0u3) */
	ul = SHL(ul) | u1;	/* (u2u1) */
	p0 = quot * v0;
	p1 = quot * v1;
	s = p0 + SHL(p1);
	ul -= s;
	borrow = (ul > ~s ? 1 : 0);
	uh -= MSH(p1) - borrow;

	if (0 != MSH(uh)) {
		/* add back v */
		quot--;
		ul += v;
		uh = 0;
	}

	/* quotient q1 */
	q1 = quot;

	/*
	 * ROUND 3:
	 * calculate q0; estimate quot = (u2u1) / v1,
	 * then set (u2u1u0) -= quot(v1v0)
	 */
	t = ul;
	quot = t / v1;
	rem = t - quot * v1;
	/* Test on v0 */
	t = SHL(rem) | u0;
	if (B == quot || (quot * v0) > t) {
		quot--;
		rem += v1;
		t = SHL(rem) | u0;
		if (rem < B && (quot * v0) > t)
			quot--;
	}

	/*
	 * multiply and subtract:
	 * (u2u1u0)" = (u2u1u0)' - quot(v1v0)
	 */
	uh = MSH(ul);			/* (0u2) */
	ul = SHL(ul) | u0;	/* (u1u0) */

	p0 = quot * v0;
	p1 = quot * v1;
	s = p0 + SHL(p1);
	ul -= s;
	borrow = (ul > ~s ? 1 : 0);
	uh -= MSH(p1) - borrow;
	if (0 != MSH(uh)) {
		/* add back v */
		quot--;
		ul += v;
		uh = 0;
	}

	/* quotient q1q0 */
	*q = SHL(q1) | LSH(quot);

	/* Remainder is in (u1u0) i.e. ul */
	*r = ul;

	/* quotient q2 (overflow) is returned */
	return q2;
}

#endif	/* HAVE_U64 */

/**
 * @brief Return greatest common divisor of two single limbs
 *
 * Returns gcd(x, y)
 *
 * Ref: Schneier 2nd ed, p245
 *
 * @param x single limb candidate #1
 * @param y single limb candidate #2
 *
 * @result return zero if x and y are zero, else gcd(x,y)
 */
limb_t sl_gcd(limb_t x, limb_t y)
{
	limb_t g;

	if (x + y == 0)
		return 0;	/* Error */

	g = y;
	while (x > 0) {
		g = x;
		x = y % x;
		y = g;
	}
	return g;
}

/**
 * @brief Test a single limb if it is a probable prime
 *
 * Returns 1 if w is a probable prime 
 * Runs t iterations of Rabin-Miller (t = 50 for DSS Standard) 
 * 
 * Uses Rabin-Miller Probabilistic Primality Test,
 * Ref: FIPS-186-2 Appendix 2.
 * Also Schneier 2nd ed p 260 & Knuth Vol 2, p 379.
 *
 * @param w Single limb probable prime to test
 * @param t Number of rounds of Rabin-Miller to apply
 *
 * @result zero if w is no prime, 1 if w is probably prime
 */
int sl_isprime(limb_t w, cnt_t t)
{

	cnt_t i, j;
	limb_t m;
	limb_t a;
	limb_t b;
	limb_t z;
	limb_t w1;
	int failed;

	/*	First check for the first t small primes */
	for (i = 0; i < SM_FIRST; i++) {
		if (w == sm_primes[i])
			return 1;	/* prime */
		if (0 == (w % sm_primes[i]))
			return 0;	/* found a divisor */
	}

	w1 = w - 1;
	/*
	 * Now do Rabin-Miller:
	 * Find a and m where w = (2^b) * m + 1.
	 * So m is odd, and 2^b is largest power of 2
	 * dividing w - 1.
	 */
	m = w1;
	b = 0;
	while (sl_iseven(m)) {
		/* Divide by 2 until m is odd */
		m >>= 1;
		b++;
	}

	/* assert((1 << b) * m + 1 == w); */

	for (i = 0; i < t; i++) {
		/* Assume fail unless passed in loop */
		failed = 1;

		/* return a random integer 1 < a <= w1 (< w) */
		a = sl_random(2, w1);

		j = 0;

		/* z = (a^m) % w */
		sl_modexp(&z, a, m, w);
		if (z == 1 || z == w1) {
			/* z = 1 or z = w - 1: passes in first test */
			failed = 0;
		} else {
			/* test up to b powers of z, too */
			do {
				j++;
				/* z = (z^2) % w */
				sl_modmul(&z, z, z, w);
				if (z == w1) {
					/* z = w - 1: passes in j-th test */
					failed = 0;
					break;
				}
				if (z == 1) {
					/* z = 1: fails */
					failed = 1;
					break;
				}
			} while (j < b);
		}

		if (failed) {
			/* w is not a prime: bail out */
			return 0;
		}
	}

	/* check for all the small primes from SM_FIRST to SM_PRIMES */
	for (i = SM_FIRST; i < SM_PRIMES; i++) {
		if (w == sm_primes[i])
			return 1;	/* is prime */
		if (0 == (w % sm_primes[i]))
			return 0;	/* found a divisor */
	}

	/* w is probably prime */
	return 1;
}

/**
 * @brief Return a single limb "RSA prime"
 *
 * Note: This isn't actually used, as single limb primes
 * would be a little too easy to guess right.
 *
 * @result A probable prime
 */
limb_t sl_rsaprime(void)
{
	limb_t p;

	p = sl_random(0, LIMBMASK) | DIGLSB;
	while (!sl_isprime(p, 50)) {
		p += 2;
	}
	return p;
}

/**
 * @brief Compute single limb exp = x^e mod m
 *
 * Computes exp = x^e mod m
 * Binary left-to-right method
 *
 * @param exp pointer to limb to receive result
 * @param x single limb x (base)
 * @param e single limb e (exponent)
 * @param m single limb m (modulus)
 *
 * @result zero on success (always!?)
 */
int sl_modexp(limb_t *exp, limb_t x, limb_t e, limb_t m)
{
	limb_t mask;
	limb_t y;	/* Temp variable */

	/* Find most significant bit in e */
	for (mask = DIGMSB; mask > 0; mask >>= 1) {
		if (e & mask)
			break;
	}

	y = x;

	for (mask >>= 1; mask > 0; mask >>= 1) {
		sl_modmul(&y, y, y, m);		/* y = (y^2) % m */
		if (e & mask)
			sl_modmul(&y, y, x, m);	/* y = (y*x) % m*/
	}

	*exp = y;
	return 0;
}

/**
 * @brief Compute single limb inverse inv = u^(-1) % v
 *
 * Computes inv = u^(-1) % v
 * Ref: Knuth Algorithm X Vol 2 p 342 
 * ignoring u2, v2, t2 and avoiding negative numbers
 *
 * @param inv pointer to limb to receive result
 * @param u single limb to inverse
 * @param v single limb modulus
 *
 * @result zero on success (always!?)
 */
int sl_modinv(limb_t *inv, limb_t u, limb_t v)
{
	limb_t u1, u3, v1, v3, t1, t3, q, w;
	int iter = 1;
	
	/* Step X1. Initialize */
	u1 = 1;
	u3 = u;
	v1 = 0;
	v3 = v;

	/* Step X2. */
	while (v3 != 0) {
		/* Step X3. */
		q = u3 / v3;	/* Divide and */
		t3 = u3 % v3;
		w = q * v1;	/* "Subtract" */
		t1 = u1 + w;
		/* Swap */
		u1 = v1;
		v1 = t1;
		u3 = v3;
		v3 = t3;
		iter = -iter;
	}

	if (iter < 0)
		*inv = v - u1;
	else
		*inv = u1;

	return 0;
}

/**
 * @brief Compute single limb a = (x * y) % mod
 *
 * Computes a = (x * y) % m
 *
 * @param a pointer to single limb to receive result
 * @param x single limb factor 1
 * @param y single limb factor 2
 * @param m single limb modulus
 *
 * @result zero on success (always!?)
 */
int sl_modmul(limb_t *a, limb_t x, limb_t y, limb_t m)
{
	static limb_t pp[2];

	/* pp[] = x * y */
	sl_mul(pp, x, y);

	/* *a = pp[] % m */
	*a = bn_mod_limb(pp, m, 2);

	/* Clean temp */
	pp[0] = pp[1] = 0;
	return 0;
}

#if	HAVE_U64
/**
 * @brief Compute double limb product of two single limbs
 *
 * Computes p[] = x * y
 * where p is two limbs (double precision) and x, y are single
 * limbs. Use double precision natively supported on this machine.
 *
 * @param p pointer to an array of two limbs receiving the result
 * @param x single limb factor #1
 * @param y single limb factor #2
 *
 * @result zero on success (always)
 */
int sl_mul(limb_t p[2], limb_t x, limb_t y)
{
	dlimb_t dd;

	dd = (dlimb_t)x * y;
	p[0] = (limb_t)dd;
	p[1] = (limb_t)(dd >> 32);
	return 0;
}

#else

/**
 * @brief Compute double limb product of two single limbs
 *
 * Computes p[] = x * y
 * Source: Arbitrary Precision Computation
 * http://numbers.computation.free.fr/Constants/constants.html
 *
 * The limbs x and y are split in halves and the four products
 * x1*y1, x0*y1, x1*y0 and x0*y0 are added shifting them to
 * their respective least significant bit position:
 * p[1] = x1*y1 + high(x0*y1 + x1*y0) + ch << 16 + cl
 * p[0] = x0*y0 + low(x0*y1 + x1*y0) << 16
 * ch = carry from adding x0*y1 + x1*y0
 * cl = carry from adding low(x0*y1 + x1*y0) << 16 to p[0]
 *
 * @param p pointer to an array of two limbs receiving the result
 * @param x single limb factor #1
 * @param y single limb factor #2
 *
 * @result zero on success (always)
 */
int sl_mul(limb_t p[2], limb_t x, limb_t y)
{
	limb_t x0, y0, x1, y1;
	limb_t t, u, carry;

	/*
	 * Split each x,y into two halves
	 *   x = x0 + B*x1
	 *   y = y0 + B*y1
	 * where B = 2^16, half the limb size
	 * Product is
	 *   xy = x0y0 + B(x0y1 + x1y0) + B^2(x1y1)
	 */
	x0 = LSH(x);
	x1 = MSH(x);
	y0 = LSH(y);
	y1 = MSH(y);

	/* Compute low part (w/o carry) */
	p[0] = x0 * y0;

	/* middle part */
	t = x0 * y1;
	u = x1 * y0;
	t += u;
	carry = (t < u ? 1 : 0);

	/*
	 * The carry will go to high half of p[1],
	 * and the high half of t will go into the
	 * into low half of p[1]
	 */
	carry = SHL(carry) + MSH(t);

	/* add low half of t to high half of p[0] */
	t = SHL(t);
	p[0] += t;
	if (p[0] < t)
		carry++;

	p[1] = x1 * y1 + carry;

	return 0;
}

#endif	/* HAVE_U64 */

/**
 * @brief Returns a single limb pseudo random
 *
 * Returns a single limb pseudo random number in the
 * range between lower and upper inclusively.
 *
 * Handles opening /dev/urandom on its own.
 *
 * @param lower single limb lower boundary for number
 * @param upper single limb upper boundary for number
 *
 * @result a pseudo random number betwee lower and upper
 */
limb_t sl_random(limb_t lower, limb_t upper)
{
	static int urandom = 0;
	limb_t r;
	limb_t mask;

	if (upper < lower) {
		fprintf(stderr, "upper < lower (%x < %x)\n",
			upper, lower);
		exit(1);
	}
	if (0 == urandom)
		urandom = open("/dev/urandom", O_RDONLY|O_BINARY);
	if (-1 == urandom) {
		fprintf(stderr, "opening /dev/urandom failed (%s)\n",
			strerror(errno));
		exit(1);
	}
	mask = DIGMSB;
	while (0 == (upper & mask))
		mask = mask | (mask >> 1);
	mask <<= 1;
	do {
		int i = read(urandom, &r, sizeof(r));
		if (-1 == i) {
			fprintf(stderr, "reading /dev/urandom failed (%s)\n",
				strerror(errno));
			exit(1);
		}
		/* discard bits that are not in the mask */
		r &= ~mask;
	} while (r < lower || r > upper);

	return r;
}

/**
 * @brief Compute division of big number by a "half digit"
 *
 * Computes q[] = u[] / v, also returns r = u[] % v
 * where q, a are big number integers of nlimb limbs each,
 * and d, r are single limbs
 *
 * Using bit-by-bit method from MSB to LSB,
 * so v must be <= HALFMASK
 * 
 * According to "Principles in PGP by Phil Zimmermann"
 *
 * @param q pointer to an array of limbs to receive the result
 * @param u pointer to an array of limbs (dividend)
 * @param v single limb (actually half limb) divisor
 * @param nlimb number of limbs in the arrays
 *
 * @result returns remainder of the division
 */
limb_t bn_div_hdig(limb_t q[], limb_t u[], limb_t v, cnt_t nlimb)
{
	limb_t mask = DIGMSB;
	limb_t r = 0;
	if (v > HALFMASK) {
		fprintf(stderr, "bn_div_hdig called with v:%x\n", v);
		exit(1);
	}

	if (0 == nlimb)
		return 0;
	if (0 == v)
		return 0;	/* Divide by zero error */

	/* Initialize quotient */
	bn_zero(q, nlimb);
	
	/* Work from MSB to LSB */
	while (nlimb > 0) {
		/* Multiply remainder by 2 */
		r <<= 1;

		/* Look at current bit */
		if (u[nlimb-1] & mask)
			r++;
		if (r >= v) {
			/* Remainder became greater than divisor */
			r -= v;
			q[nlimb-1] |= mask;
		}

		/* next bit */
		mask >>= 1;
		if (0 != mask)
			continue;

		/* next limb */
		--nlimb;
		mask = DIGMSB;
	}
	return r;
}

/**
 * @brief Compute single limb remainder of bignum % single limb
 *
 * Computes r = u[] % v
 * where a is a big number integer of nlimb
 * and r, v are single limbs, using bit-by-bit
 * method from MSB to LSB.
 * 
 * Ref:
 *   Derived from principles in PGP by Phil Zimmermann
 * Note:
 *   This method will only work until r <<= 1 overflows.
 *   i.e. for d < DIGMSB, but we keep HALF_DIGIT
 *   limit for safety, and also because we don't
 *   have a 32nd bit.
 *
 * @param u pointer to big number to divide
 * @param v single limb (actually half limb) modulus
 * @param nlimb number of limbs in the array
 *
 * @result returns remainder of the division
 */
limb_t bn_mod_hdig(limb_t u[], limb_t v, cnt_t nlimb)
{
	limb_t mask;
	limb_t r;

	if (0 == nlimb)
		return 0;
	if (0 == v)
		return 0;	/* Divide by zero error */

	if (v > HALFMASK) {
		fprintf(stderr, "bn_mod_hdig called with v:%x\n", v);
		exit(1);
	}

	/* Work from left to right */
	mask = DIGMSB;
	r = 0;
	while (nlimb > 0) {
		/* Multiply remainder by 2 */
		r <<= 1;

		/* Look at current bit */
		if (u[nlimb-1] & mask)
			r++;

		if (r >= v)
			/* Remainder became greater than divisor */
			r -= v;

		/* next bit */
		mask >>= 1;
		if (0 != mask)
			continue;

		/* next limb */
		--nlimb;
		mask = DIGMSB;
	}
	return r;
}

/**
 * @brief Addition of two bignum arrays
 *
 * Computes w[] = u[] + v[]
 * where w, u, v are big number integers of nlimb limbs each.
 * Returns carry, i.e. w[nlimb], as 0 or 1.
 * 
 * Ref: Knuth Vol 2 Ch 4.3.1 p 266 Algorithm A.
 *
 * @param w pointer to array of limbs to receive the result
 * @param u pointer to array of limbs (addend #1)
 * @param v pointer to array of limbs (addend #2)
 * @param nlimb number of limbs in the arrays
 *
 * @result returns the carry, i.e. w[nlimb], as 0 or 1
 */
limb_t bn_add(limb_t w[], limb_t u[], limb_t v[], cnt_t nlimb)
{
	limb_t carry;
	cnt_t j;

	for (j = 0, carry = 0; j < nlimb; j++) {
		/*
		 * add limbs w[j] = u[j] + v[j] + carry;
		 * set carry = 1 if carry (overflow) occurs
		 */
		w[j] = u[j] + carry;
		carry = (w[j] < carry ? 1 : 0);

		w[j] = w[j] + v[j];
		if (w[j] < v[j])
			carry++;
	}

	/* w[n] = carry */
	return carry;
}

/**
 * @brief Subtraction of two bignum arrays
 *
 * Calculates w[] = u[] - v[] where u[] >= v[]
 * w, u, v are big number integers of nlimb limbs each
 * Returns 0 if ok, or 1 if v was greater than u.
 * 
 * Ref: Knuth Vol 2 Ch 4.3.1 p 267 Algorithm S.
 *
 * @param w pointer to array of limbs to receive the result
 * @param u pointer to array of limbs (minuend)
 * @param v pointer to array of limbs (subtrahend)
 * @param nlimb number of limbs in the arrays
 *
 * @result zero on success, 1 if v was greater than u
 */
limb_t bn_sub(limb_t w[], limb_t u[], limb_t v[], cnt_t nlimb)
{
	limb_t borrow;
	cnt_t j;

	for (j = 0, borrow = 0; j < nlimb; j++) {
		/*
		 * Subtract limbs w[j] = u[j] - v[j] - borrow;
		 * set borrow = 1 if borrow occurs
		 */
		w[j] = u[j] - borrow;
		borrow = (w[j] > ~borrow ? 1 : 0);

		w[j] = w[j] - v[j];
		if (w[j] > ~v[j])
			borrow++;
	}

	/* borrow should be 0, if u >= v */
	return borrow;
}

/**
 * @brief Product of two bignum arrays
 *
 * Computes product w[] = u[] * v[]
 * where u, v are big number integers of nlimb each
 * and w is a big number integer of 2*nlimb limbs.
 * 
 * Ref: Knuth Vol 2 Ch 4.3.1 p 268 Algorithm M.
 *
 * @param w pointer to array of limbs to receive the result
 * @param u pointer to array of limbs (factor #1)
 * @param v pointer to array of limbs (factor #2)
 * @param nlimb number of limbs in the arrays
 *
 * @result zero on success (always!?)
 */
int bn_mul(limb_t w[], limb_t u[], limb_t v[], cnt_t nlimb)
{
	limb_t t[2];
	limb_t carry;
	cnt_t i, j, m, n;

	m = n = nlimb;

	/* zero result */
	bn_zero(w, 2*nlimb);

	for (j = 0; j < n; j++) {
		/* zero multiplier? */
		if (0 == v[j]) {
			w[j+m] = 0;
			continue;
		}
		/* Initialize i */
		carry = 0;
		for (i = 0; i < m; i++) {
			/*
			 * Multiply and add:
			 * t = u[i] * v[j] + w[i+j] + carry
			 */
			sl_mul(t, u[i], v[j]);

			t[0] += carry;
			if (t[0] < carry)
				t[1]++;
			t[0] += w[i+j];
			if (t[0] < w[i+j])
				t[1]++;

			w[i+j] = t[0];
			carry = t[1];
		}	
		w[j+m] = carry;
	}

	return 0;
}

/**
 * @brief Shift left a bignum by a number of bits (less than LIMBBITS)
 *
 * Computes a[] = b[] << x
 * Where a and b are big number integers of nlimb each.
 * The shift count must be less than LIMBBITS
 *
 * @param a pointer to array of limbs to receive the result
 * @param b pointer to array of limbs to shift left
 * @param x number of bits to shift (must be less than LIMBBITS)
 * @param nlimb number of limbs in the arrays
 *
 * @result returns a single limb "carry", i.e. bits that came out left
 */
limb_t bn_shl(limb_t a[], limb_t b[], cnt_t x, cnt_t nlimb)
{
	cnt_t i, y;
	limb_t carry, temp;

	if (0 == nlimb)
		return 0;

	if (0 == x) {
		/* no shift at all */
		if (a != b)
			bn_copy(a, b, nlimb);
		return 0;
	}

	/* check shift amount */
	if (x >= LIMBBITS) {
		fprintf(stderr, "bn_shl() called with x >= %d\n", LIMBBITS);
		exit(1);
		return 0;
	}

	y = LIMBBITS - x;
	carry = 0;
	for (i = 0; i < nlimb; i++) {
		temp = b[i] >> y;
		a[i] = (b[i] << x) | carry;
		carry = temp;
	}

	return carry;
}

/**
 * @brief Shift right a bignum by a number of bits (less than LIMBBITS)
 *
 * Computes a[] = b[] >> x
 * Where a and b are big number integers of nlimb each.
 * The shift count must be less than LIMBBITS
 *
 * @param a pointer to array of limbs to receive the result
 * @param b pointer to array of limbs to shift right
 * @param x number of bits to shift (must be less than LIMBBITS)
 * @param nlimb number of limbs in the arrays
 *
 * @result returns a single limb "carry", i.e. bits that came out right
 */
limb_t bn_shr(limb_t a[], limb_t b[], cnt_t x, cnt_t nlimb)
{
	cnt_t i, y;
	limb_t carry, temp;

	if (0 == nlimb)
		return 0;

	if (0 == x) {
		/* no shift at all */
		if (a != b)
			bn_copy(a, b, nlimb);
		return 0;
	}

	/* check shift amount */
	if (x >= LIMBBITS) {
		fprintf(stderr, "bn_shr() called with x >= %d\n", LIMBBITS);
		exit(1);
	}

	y = LIMBBITS - x;
	carry = 0;
	i = nlimb;
	while (i-- > 0) {
		temp = b[i] << y;
		a[i] = (b[i] >> x) | carry;
		carry = temp;
	}

	return carry;
}

/**
 * @brief Check a quotient for overflow
 *
 * Returns 1 if quot is too big,
 * i.e. if (quot * Vn-2) > (b.rem + Uj+n-2)
 * Returns 0 if ok
 *
 * @param quot quotient under test
 * @param rem remainder
 * @param
 *
 * @result zero on success
 */
static int quot_overflow(limb_t quot, limb_t rem, limb_t v, limb_t u)
{
	limb_t t[2];

	sl_mul(t, quot, v);
	if (t[1] < rem)
		return 0;
	if (t[1] > rem)
		return 1;
	if (t[0] > u)
		return 1;

	return 0;
}

/**
 * @brief Compute quotient and remainder of bignum division
 *
 * Computes quotient q[] = u[] / v[]
 * and remainder r[] = u[] % v[]
 * where q, r, u are big number integers of ulimb limbs,
 * and the divisor v of vlimb limbs.
 * 
 * Ref: Knuth Vol 2 Ch 4.3.1 p 272 Algorithm D.
 *
 * @param q pointer to array of limbs to receive quotient
 * @param r pointer to array of limbs to receive remainder
 * @param u pointer to array of limbs (dividend)
 * @param ulimb number of limbs in the q, r, u arrays
 * @param v pointer to array of limbs (divisor)
 * @param vlimb number of limbs in the v array
 *
 * @result zero on success, LIMBASK on division by zero
 */
int bn_div(limb_t q[], limb_t r[], limb_t u[], limb_t v[],
	cnt_t ulimb, cnt_t vlimb)
{
	static limb_t qq[BN_MAXSIZE];
	static limb_t uu[BN_MAXSIZE];
	static limb_t vv[BN_MAXSIZE];
	limb_t mask;
	limb_t overflow;
	limb_t quot;
	limb_t rem;
	limb_t t[2];
	limb_t *ww;
	cnt_t n, m, i, j, shift;
	int ok, cmp;

	/* find size of v */
	n = bn_sizeof(v, vlimb);

	/* Catch special cases */
	if (0 == n)
		return LIMBMASK;	/* Error: divide by zero */

	if (1 == n) {
		/* Use short division instead */
		r[0] = bn_div_limb(q, u, v[0], ulimb);
		return 0;
	}

	/* find size of u */
	m = bn_sizeof(u, ulimb);

	if (m < n) {
		/* v > u: just set q = 0 and r = u */
		bn_zero(q, ulimb);
		bn_copy(r, u, ulimb);
		return 0;
	}

	if (m == n) {
		/* u and v are the same length: compare them */
		cmp = bn_cmp(u, v, (unsigned int)n);
		if (0 == cmp) {
			/* v == u: set q = 1 and r = 0 */
			bn_limb(q, 1, ulimb);
			bn_zero(r, ulimb);
			return 0;
		}
		if (cmp < 0) {
			/* v > u: set q = 0 and r = u */
			bn_zero(q, ulimb);
			bn_copy(r, u, ulimb);
			return 0;
		}
	}

	/* m greater than or equal to n */
	m -= n;

	/* clear quotient qq */
	bn_zero(qq, ulimb);

	/*
	 * Normalize v: requires high bit of v[n-1] to be set,
	 * so find most significant bit, then shift left
	 */
	mask = DIGMSB;
	for (shift = 0; shift < LIMBBITS; shift++) {
		if (v[n-1] & mask)
			break;
		mask >>= 1;
	}

	/* normalize vv from v */
	overflow = bn_shl(vv, v, shift, n);

	/* copy normalized dividend u into remainder uu */
	overflow = bn_shl(uu, u, shift, n + m);

	/* new limb u[m+n] */
	t[0] = overflow;

	j = m + 1;
	while (j-- > 0) {
		/* quot = (b * u[j+n] + u[j+n-1]) / v[n-1] */
		ok = 0;

		/* This is Uj+n */
		t[1] = t[0];
		t[0] = uu[j+n-1];

		overflow = sl_div(&quot, &rem, t, vv[n-1]);

		if (overflow) {
			/* quot = b */
			quot = LIMBMASK;
			rem = uu[j+n-1] + vv[n-1];
			if (rem < vv[n-1])
				ok = 1;
		}
		if (0 == ok && quot_overflow(quot, rem, vv[n-2], uu[j+n-2])) {
			/* quot * v[n-2] > b * rem + u[j+n-2] */
			quot--;
			rem += vv[n-1];
			if (rem >= vv[n-1])
				if (quot_overflow(quot, rem, vv[n-2], uu[j+n-2]))
					quot--;
		}

		/* multiply and subtract vv[] * quot */
		ww = &uu[j];

		if (0 == quot) {
			overflow = 0;
		} else {
			/* quot is non zero */
			limb_t tt[2];
			limb_t borrow;

			for (i = 0, borrow = 0; i < n; i++) {
				sl_mul(tt, quot, vv[i]);
				ww[i] -= borrow;
				borrow = (ww[i] > ~borrow ? 1 : 0);

				ww[i] -= tt[0];
				if (ww[i] > ~tt[0])
					borrow++;
				borrow += tt[1];
			}

			/*
			 * w[n] is not in array w[0..n-1]:
			 * subtract final borrow
			 */
			overflow = t[1] - borrow;
		}

		/* test for remainder */
		if (overflow) {
			quot--;
			/* add back if mul/sub was negative */
			overflow = bn_add(ww, ww, vv, n);
		}

		qq[j] = quot;

		/* u[j+n] for next round */
		t[0] = uu[j+n-1];
	}

	/* clear uu[] limbs from n to n+m */
	for (j = n; j < m+n; j++)
		uu[j] = 0;

	/* denormalize remainder */
	bn_shr(r, uu, shift, n);

	/* copy quotient */
	bn_copy(q, qq, n + m);

	/* clear temps */
	bn_zero(qq, n);
	bn_zero(uu, n);
	bn_zero(vv, n);
	return 0;
}

/**
 * @brief Compute remainder of bignum division (modulo)
 *
 * Calculates r[] = u[] % v[]
 * where r, v are big number integers of length vlimb
 * and u is a big number integer of length ulimb.
 * r may overlap v.
 * 
 * Note that r here is only vlimb long, 
 * whereas in bn_div it is ulimb long.
 * 
 * Use remainder from bn_div function.
 *
 * @param r pointer to array of limbs to receive remainder
 * @param u pointer to array of limbs (dividend)
 * @param ulimb number of limbs in the u array
 * @param v pointer to array of limbs (divisor)
 * @param vlimb number of limbs in the r and v array
 *
 * @result zero on success, LIMBASK on division by zero
 */
limb_t bn_mod(limb_t r[], limb_t u[], cnt_t ulimb, limb_t v[], cnt_t vlimb)
{
	static limb_t qq[2*BN_MAXSIZE];
	static limb_t rr[2*BN_MAXSIZE];
	limb_t d0;

	/* rr[] = u[] % v[n] */
	d0 = bn_div(qq, rr, u, v, ulimb, vlimb);

	/* copy vlimb limbs of remainder */
	bn_copy(r, rr, vlimb);

	/* zero temps */
	bn_zero(rr, ulimb);
	bn_zero(qq, ulimb);

	return d0;
}

/**
 * @brief Compute greatest common divisor
 *
 * Computes g = gcd(x, y)
 * Reference: Schneier
 *
 * @param g pointer to array of limbs to receive the gcd
 * @param x pointer to array of limbs (candidate #1)
 * @param y pointer to array of limbs (candidate #2)
 * @param nlimb number of limbs in the arrays
 *
 * @result zero on succes (always)
 */
int bn_gcd(limb_t g[], limb_t x[], limb_t y[], cnt_t nlimb)
{	
	static limb_t yy[BN_MAXSIZE];
	static limb_t xx[BN_MAXSIZE];
	
	bn_copy(xx, x, nlimb);
	bn_copy(yy, y, nlimb);

	/* g = y */
	bn_copy(g, yy, nlimb);
	
	/* while (x > 0) { */
	while (0 != bn_cmp_limb(xx, 0, nlimb)) {
		/* g = x */
		bn_copy(g, xx, nlimb);
		/* x = y % x */
		bn_mod(xx, yy, nlimb, xx, nlimb);
		/* y = g */
		bn_copy(yy, g, nlimb);
	}

	bn_zero(xx, nlimb);
	bn_zero(yy, nlimb);
	
	/* gcd is left in g */
	return 0;
}

/**
 * @brief Compute modular exponentiation of bignums
 *
 * Computes y[] = (x[]^e[]) % m[]
 * Binary MSB to LSB method
 *
 * @param y pointer to array of limbs to receive the result
 * @param x pointer to array of limbs (base)
 * @param e pointer to array of limbs (exponent)
 * @param m pointer to array of limbs (modulus)
 * @param nlimb number of limbs in the arrays
 *
 * @result zero on success, -1 on error (nlimb is zero)
 */
int bn_modexp(limb_t y[], limb_t x[], limb_t e[], limb_t m[], cnt_t nlimb)
{
	limb_t mask;
	cnt_t n;
	
	if (nlimb == 0)
		return -1;

	/* Find second-most significant bit in e */
	n = bn_sizeof(e, nlimb);
	for (mask = DIGMSB; 0 != mask; mask >>= 1) {
		if (e[n-1] & mask)
			break;
	}
	/* next bit, because we start off with y[] == x[] */
	mask >>= 1;
	if (0 == mask) {
		mask = DIGMSB;
		n--;
	}

	/* y[] = x[] */
	bn_copy(y, x, nlimb);

	while (n > 0) {
		/* y[] = (y[] ^ 2) % m[] */
		bn_modmul(y, y, y, m, nlimb);

		if (e[n-1] & mask)
			/* y[] = (y[] * x[]) % m[] */
			bn_modmul(y, y, x, m, nlimb);

		/* next bit */
		mask >>= 1;
		if (0 == mask) {
			mask = DIGMSB;
			n--;
		}
	}

	return 0;
}

/**
 * @brief Compute modular product of two bignums
 *
 * Computes a[] = (x[] * y[]) % m[]
 * where a, x, y and m are big numbers of nlimb length
 *
 * @param a pointer to array of limbs to receive the result
 * @param x pointer to array of limbs (factor #1)
 * @param y pointer to array of limbs (factor #2)
 * @param m pointer to array of limbs (modulus)
 * @param nlimb number of limbs in the arrays
 *
 * @result zero on success, LIMBMASK if m was zero (division by zero)
 */
limb_t bn_modmul(limb_t a[], limb_t x[], limb_t y[], limb_t m[], cnt_t nlimb)
{
	static limb_t pp[2*BN_MAXSIZE];
	limb_t d0;

	/* pp[] = x[] * y[] (NB: double size pp[]) */
	bn_mul(pp, x, y, nlimb);

	/* a[] = pp[] % m[] */
	d0 = bn_mod(a, pp, 2*nlimb, m, nlimb);

	/* zero temp */
	bn_zero(pp, 2*nlimb);

	return d0;
}

/**
 * @brief Compute modular inverse
 *
 * Computes inv[] = u[]^(-1) % v[]
 * Ref: Knuth Algorithm X Vol 2 p 342
 * ignoring u2, v2, t2 and avoiding negative numbers.
 *
 * @param inv pointer to array of limbs receiving the result
 * @param u pointer to array of limbs (candidate)
 * @param v pointer to array of limbs (modulus)
 * @param nlimb number of limbs in the arrays
 *
 * @result zero on success
 */
int bn_modinv(limb_t inv[], limb_t u[], limb_t v[], cnt_t nlimb)
{
	/* Allocate temp variables */
	static limb_t u1[BN_MAXSIZE];
	static limb_t u3[BN_MAXSIZE];
	static limb_t v1[BN_MAXSIZE];
	static limb_t v3[BN_MAXSIZE];
	static limb_t t1[BN_MAXSIZE];
	static limb_t t3[BN_MAXSIZE];
	static limb_t q[BN_MAXSIZE];
	static limb_t w[2*BN_MAXSIZE];
	int iter;

	/* Step X1. Initialize */
	bn_limb(u1, 1, nlimb);	/* u1 = 1 */
	bn_limb(v1, 0, nlimb);	/* v1 = 0 */
	bn_copy(u3, u, nlimb);	/* u3 = u */
	bn_copy(v3, v, nlimb);	/* v3 = v */

	/* remember odd/even iterations */
	iter = 1;

	/* Step X2. Loop while v3 != 0 */
	while (0 != bn_cmp_limb(v3, 0, nlimb)) {
		/* Step X3. Divide and "Subtract" */
		/* q = u3 / v3, t3 = u3 % v3 */
		bn_div(q, t3, u3, v3, nlimb, nlimb);
		/* w = q * v1 */
		bn_mul(w, q, v1, nlimb);
		/* t1 = u1 + w */
		bn_add(t1, u1, w, nlimb);

		/* Swap u1 <= v1 <= t1 */
		bn_copy(u1, v1, nlimb);
		bn_copy(v1, t1, nlimb);

		/* Swap u3 <= v3 <= t3 */
		bn_copy(u3, v3, nlimb);
		bn_copy(v3, t3, nlimb);

		iter ^= 1;
	}

	if (iter)
		bn_copy(inv, u1, nlimb);	/* inv = u1 */
	else
		bn_sub(inv, v, u1, nlimb);	/* inv = v - u1 */

	/* clear temp vars */
	bn_zero(u1, nlimb);
	bn_zero(v1, nlimb);
	bn_zero(t1, nlimb);
	bn_zero(u3, nlimb);
	bn_zero(v3, nlimb);
	bn_zero(t3, nlimb);
	bn_zero(q, nlimb);
	bn_zero(w, 2*nlimb);

	return 0;
}

/**
 * @brief Compute square root (and fraction) of a bignum
 *
 * Compute q[] = sqrt(u[]),
 * where q and u are big number integers of nlimb limbs
 *
 * Method according to sqrt.html of 2001-08-15:
 * Act on bytes from MSB to LSB, counting the number of times
 * that we can subtract consecutive odd numbers starting with
 * 1, 3, 5. Just uses add, subtract, shift and comparisons.
 *
 * The pointer r can be NULL if caller is not interested in
 * the (partial) fraction.
 *
 * @param q pointer to array of limbs to receive the result (integer)
 * @param r pointer to array of limbs to receive the result (fraction)
 * @param u pointer to array of limbs (square)
 * @param rlimb number of limbs in the q and r arrays
 * @param ulimb number of limbs in the u array
 *
 * @result zero on success
 */
int bn_sqrt(limb_t q[], limb_t r[], limb_t u[], cnt_t rlimb, cnt_t ulimb)
{
	static limb_t step[BN_MAXSIZE];
	static limb_t accu[BN_MAXSIZE];
	static limb_t w[2*BN_MAXSIZE];
	limb_t d;
	cnt_t m, n;
	cnt_t shift;

	bn_zero(q, ulimb);
	bn_limb(step, 1, BN_MAXSIZE);
	bn_limb(accu, 0, BN_MAXSIZE);
	n = bn_sizeof(u, ulimb);

	/* determine first non-zero byte from MSB to LSB */
	if (0 != (u[n-1] >> 24)) {
		shift = 32;
	} else if (0 != (u[n-1] >> 16)) {
		shift = 24;
	} else if (0 != (u[n-1] >> 8)) {
		shift = 16;
	} else {
		shift = 8;
	}

	m = 1;
	while (n-- > 0) {
		while (shift > 0) {
			/* shift accu one byte left */
			bn_shl(accu, accu, 8, m+1);

			/* shift for next byte from u[] */
			shift -= 8;
			accu[0] |= (u[n] >> shift) & 0xff;

			/* digit = 0 */
			d = 0;
			/* subtract consecutive odd numbers step[] until overflow */
			for (d = 0; bn_cmp(step, accu, m+1) <= 0; d++) {
				bn_sub(accu, accu, step, m+1);
				bn_add_limb(step, step, 2, m+1);
			}

			/* put digit into result */
			bn_shl(q, q, 4, m);
			q[0] |= d;

			/* step[] = 2 * q[] * 16 + 1 */
			bn_shl(step, q, 5, m+1);
			bn_add_limb(step, step, 1, m+1);
		}
		shift = 32;
		if (0 == (n & 1))
			m++;
	}

	/* Caller does not want to know the fraction? */
	if (NULL == r)
		return 0;

	/* nothing left to do if remainder is zero */
	if (0 == bn_cmp_limb(accu, 0, ulimb)) {
		bn_zero(r, rlimb);
		return 0;
	}

	/* Start off with the integer part */
	bn_zero(w, 2*BN_MAXSIZE);
	bn_copy(w, q, ulimb);

	n = rlimb * (LIMBBITS / 4);
	while (n-- > 0) {
		/* shift accu one byte left */
		bn_shl(accu, accu, 8, rlimb);

		/* subtract consecutive odd numbers step[] until overflow */
		for (d = 0; bn_cmp(step, accu, rlimb) <= 0; d++) {
			bn_sub(accu, accu, step, rlimb);
			bn_add_limb(step, step, 2, rlimb);
		}

		/* put digit into result */
		bn_shl(w, w, 4, rlimb);
		w[0] |= d;

		/* step[] = 2 * w[] * 16 + 1 */
		bn_shl(step, w, 5, rlimb);
		bn_add_limb(step, step, 1, rlimb);
	}

	/* copy remainder */
	bn_copy(r, w, rlimb);
	return 0;
}

/**
 * @brief Test a bignum for being a probable prime
 *
 * Returns 1 if p is a probable prime , else 0.
 * Carries out t iterations (use t = 50 for DSS Standard)
 *
 * Uses Rabin-Miller Probabilistic Primality Test,
 * Ref: FIPS-186-2 Appendix 2.
 * Also Schneier 2nd ed p 260 & Knuth Vol 2, p 379.
 *
 * @param p pointer to an array of limbs (candidate)
 * @param nlimb number of limbs in the array
 * @param t number of tries for Rabin-Miller
 *
 * @result zero if p is not prime, 1 if p is probably prime
 */
int bn_isprime(limb_t p[], cnt_t nlimb, cnt_t t)
{
	static limb_t m[BN_MAXSIZE];
	static limb_t a[BN_MAXSIZE];
	static limb_t z[BN_MAXSIZE];
	static limb_t p1[BN_MAXSIZE];
	cnt_t b, i, j;
	int failed, isprime;

	/* check the obvious */
	if (bn_iseven(p, nlimb))
		return 0;

	/* First check for the first SM_FIRST small primes */
	for (i = 0; i < SM_FIRST; i++) {
		if (0 == bn_cmp_limb(p, sm_primes[i], nlimb))
			return 1;	/* Found prime */
		if (0 == bn_mod_hdig(p, sm_primes[i], nlimb))
			return 0;	/* Failed */
	}

	/* Keep a copy of p1 = p - 1 for comparisons */
	bn_sub_limb(p1, p, 1, nlimb);

	/*
	 * Now do Rabin-Miller:
	 * Find a and m where p = (2 ^ b) * m + 1,
	 * m is odd and 2^b is largest power of 2 dividing p - 1
	 */

	/* m = p - 1 */
	bn_copy(m, p1, nlimb);

	/* number of shifts b = 0 */
	b = 0;
	while (bn_iseven(m, nlimb)) {
		/* divide by 2 until m is odd */
		bn_shr(m, m, 1, nlimb);
		/* and count in b */
		b++;
	}

	isprime = 1;
	for (i = 0; i < t; i++) {
		/* assume failure */
		failed = 1;
		
		/* make a random integer a: 1 < a <= p1 (a < p) */
		bn_zero(a, nlimb);
		do {
			a[0] = sl_random(2, LIMBMASK);
		} while (bn_cmp(a, p1, nlimb) > 0);

		/* z[] = (a[] ^ m[]) % p[] */
		bn_modexp(z, a, m, p, nlimb);

		j = 0;
		if (0 == bn_cmp_limb(z, 1, nlimb)) {
			/* z == 1 : pass */
			failed = 0;
		} else if (0 == bn_cmp(z, p1, nlimb)) {
			/* z == p - 1 : pass */
			failed = 0;
		} else {
			do {
				/* z = (z^2) % p */
				bn_modmul(z, z, z, p, nlimb);
				if (0 == bn_cmp(z, p1, nlimb)) {
					/* z = p - 1 : pass */
					failed = 0;
					break;
				}
				if (0 == bn_cmp_limb(z, 1, nlimb)) {
					/* z = 1 : failed */
					failed = 1;
					break;
				}
			} while (j++ < b);
		}

		if (failed) {
			/* p is not a prime : break out */
			isprime = 0;
			break;
		}
		/* p is still a possible prime: continue tests */
	}

	/* check for all the small primes SM_FIRST to SM_PRIMES */
	for (i = SM_FIRST; i < SM_PRIMES; i++) {
		if (0 == bn_cmp_limb(p, sm_primes[i], nlimb))
			break;	/* Found prime */
		if (0 == bn_mod_hdig(p, sm_primes[i], nlimb))
			isprime = 0;	/* found a factor! */
	}

	/* clean temp vars */
	bn_zero(m, nlimb);
	bn_zero(a, nlimb);
	bn_zero(z, nlimb);
	bn_zero(p1, nlimb);

	return isprime;
}

/**
 * @brief Return a pseudo random bignum
 *
 * Fills array r[] with nlimb pseudo random limbs
 *
 * @param r pointer to array of limbs to receive result
 * @param nlimb number of limbs in the array
 *
 * @result zero on success (always)
 */
int bn_random(limb_t r[], cnt_t nlimb)
{
	while (nlimb-- > 0)
		r[nlimb] = sl_random(0, LIMBMASK);
	return 0;
}

/**
 * @brief Generate a pseudo random prime number for RSA
 *
 * Generates a prime number for the RSA algorithm (i.e. p, q)
 *
 * The least and most significant bits are set to 1,
 * in order to save time testing for primality (avoid shifts).
 *
 * @param p pointer to an array of limbs to receive (probable) prime
 * @param nlimb number of limbs in array
 * @param t number of tries for the Rabin-Miller test
 *
 * @result zero on success (always)
 */
int bn_rsaprime(limb_t p[], cnt_t nlimb, cnt_t t)
{
	cnt_t i;

	do {
		bn_zero(p, BN_MAXSIZE);
		for (i = 0; i < nlimb; i++)
			p[i] = sl_random(0, LIMBMASK);
		p[0] |= DIGLSB;
		p[nlimb-1] |= DIGMSB + (DIGMSB >> 1);

		while (!bn_isprime(p, nlimb, t)) {
			limb_t a = sl_random(0, LIMBMASK) * 2;
			bn_add_limb(p, p, a, nlimb);
			if (DIGMSB2 != (p[nlimb-1] & DIGMSB2))
				break;
		}
	} while (DIGMSB2 != (p[nlimb-1] & DIGMSB2));
	DBG(("\n"));

	return 0;
}

/**
 * @brief Generate a RSA key pair
 *
 * Generate an RSA key pair (n[], e[]), (d[])
 * where the key size is RSALIMBS limbs (hardcoded)
 *
 * @param n pointer to an array of RSALIMBS limbs (public key)
 * @param e pointer to an array of RSALIMBS limbs (exponent)
 * @param d pointer to an array of RSALIMBS limbs (private key)
 * @param f which fermat prime to use as exponent (0 to 4)
 *
 * @result zero on success
 */
int bn_rsakey(limb_t n[], limb_t e[], limb_t d[], cnt_t f)
{
	static const limb_t fermat_primes[] = { 3, 5, 17, 257, 65537 };
	static limb_t p[RSALIMBS];
	static limb_t q[RSALIMBS];
	static limb_t g[RSALIMBS];
	static limb_t h[RSALIMBS];
	static limb_t L[2*RSALIMBS];
	cnt_t nlimb;

	/*
	 * Set e[] to the chosen Fermat number
	 * Why do we use fermat primes anyway?
	 */
	if (f > 4)
		f = 4;
	bn_zero(e, RSALIMBS);
	e[0] = fermat_primes[f];

	/* Generate a random prime p somewhere between
	 * RSALIMBS_MIN and approx. half length of modulus
	 */
	info("/p");
	bn_zero(p, RSALIMBS);
	nlimb = sl_random(RSALIMBS_MIN, RSALIMBS_MAX);
	do {
		bn_rsaprime(p, nlimb, 50);
		/* verify that gcd(p-1,e) is 1 */
		bn_sub_limb(L, p, 1, nlimb);
		bn_gcd(g, L, e, nlimb);
	} while (0 != bn_cmp_limb(g, 1, nlimb));
	DBG(("p: 0x%s\n", bn2x(p, nlimb)));

	/* generate a random prime q to balance modulus length */
	info("/q");
	bn_zero(q, RSALIMBS);
	nlimb = RSALIMBS - nlimb;
	do {
		bn_rsaprime(q, nlimb, 50);
		/* verify that gcd(q-1,e) is 1 */
		bn_sub_limb(L, q, 1, nlimb);
		bn_gcd(g, L, e, nlimb);
	} while (0 != bn_cmp_limb(g, 1, nlimb));
	DBG(("q: 0x%s\n", bn2x(q, nlimb)));

	/*
	 * Compute the modulus n = p * q
	 * use L as a double size temp
	 */

	/* L = p * q */
	bn_mul(L, p, q, RSALIMBS);

	/* copy RSALIMBS limbs */
	bn_copy(n, L, RSALIMBS);
	DBG(("n: 0x%s\n", bn2x(n, RSALIMBS)));

	/* L = (p - 1) * (q - 1) */
	bn_sub_limb(g, p, 1, RSALIMBS);
	bn_sub_limb(h, q, 1, RSALIMBS);
	bn_mul(L, g, h, RSALIMBS);
	DBG(("L: 0x%s\n", bn2x(L, RSALIMBS)));

	/* Compute private key d = (e^-1) % L */
	bn_modinv(d, e, L, RSALIMBS);
	DBG(("d: 0x%s\n", bn2x(d, RSALIMBS)));
	DBG(("e: 0x%s\n", bn2x(e, 1)));

	/* clear temp vars */
	bn_zero(p, RSALIMBS);
	bn_zero(q, RSALIMBS);
	bn_zero(g, RSALIMBS);
	bn_zero(h, RSALIMBS);
	bn_zero(L, 2*RSALIMBS);

	info(" - ");
	return 0;
}

typedef struct rsa_s {
	uint32_t bits;
	limb_t n[BN_MAXSIZE];
	limb_t e[BN_MAXSIZE];
	limb_t d[BN_MAXSIZE];
}	rsa_t;

static rsa_t mykey;

#if	!defined(BN_TEST) || (BN_TEST == 0)

/**
 * @brief Setup our RSA keypair; either load or create it
 *
 * Try to load the public key, modulus and private key from
 * the file store/.mykey, where we keep it in machine native
 * format (limbs, i.e. 32 bit values).
 * If the file isn't there, create a new keypair and write
 * it to this file.
 *
 * @result zero on success
 */
int rsa_setup(void)
{
	char filename[MAXPATHLEN];
	int fd;

	mkdir(g_conf->storepath, 0700);
	bn_snprintf(filename, sizeof(filename),
		"%s/.mykey", g_conf->storepath);
	memset(&mykey, 0, sizeof(mykey));
	/* hardcoded for now */
	mykey.bits = RSABITS;
	fd = open(filename, O_RDONLY);
	if (-1 == fd) {
		info("creating %u bits RSA keypair", RSABITS);
		bn_rsakey(mykey.n, mykey.e, mykey.d, 5);
		fd = open(filename, O_CREAT|O_TRUNC|O_RDWR|O_BINARY, 0600);
		if (-1 == fd) {
			return -1;
		}
		write(fd, mykey.n, RSALIMBS * sizeof(limb_t));
		write(fd, mykey.e, RSALIMBS * sizeof(limb_t));
		write(fd, mykey.d, RSALIMBS * sizeof(limb_t));
		close(fd);
	} else {
		read(fd, mykey.n, RSALIMBS * sizeof(limb_t));
		read(fd, mykey.e, RSALIMBS * sizeof(limb_t));
		read(fd, mykey.d, RSALIMBS * sizeof(limb_t));
		close(fd);
	}

	return 0;
}
#endif	/* !defined(BN_TEST) */

/**
 * @brief Allocate and initialize a RSA context
 *
 * @param rsa pointer to a pointer to receive RSA context
 *
 * @result zero on success, -1 on failure (errno = ENOMEM)
 */
int rsa_init(void **rsa)
{
	*rsa = calloc(1, sizeof(rsa_t));
	if (NULL != *rsa) {
		memcpy(*rsa, &mykey, sizeof(rsa_t));
		return 0;
	}
	return -1;
}

/**
 * @brief Destroy and free a RSA context
 *
 * @param rsa pointer to a pointer that was allocated by rsa_init
 *
 * @result zero on success, -1 and errno = EINVAL if called with NULL
 */
int rsa_destroy(void **rsa)
{
	if (NULL == rsa) {
		errno = EINVAL;
		return -1;
	}
	if (NULL != *rsa) {
		memset(rsa, 0, RSALIMBS * sizeof(limb_t));
		free(*rsa);
	}
	*rsa = NULL;
	return 0;
}

/**
 * @brief Extract the public key from a RSA context
 *
 * Extract the public key from a RSA context and copy it to the
 * specified buffer.
 * If this function is called with either of rsa or pubkey NULL,
 * or with pubkey_size 0, the returned value is the number of
 * bytes required to hold a copy of the public key.
 * Otherwise the return value is zero, if the buffer can
 * receive the required number of bytes for the entire public key.
 * The return value is -1 and errno set to EINVAL if the buffer
 * is too small to receive the public key.
 *
 * @param rsa pointer to a RSA context
 * @param pubkey pointer to a buffer to receive the public key
 * @param pubkey_size number of bytes available in the buffer
 *
 * @result zero on success, size of public key when called with NULL, -1 error
 */
int rsa_get_public_key(void *rsa, void *pubkey, int pubkey_size)
{
	rsa_t *r = (rsa_t *)rsa;
	uint32_t *pkey = (uint32_t *)pubkey;
	int size = RSABITS / 8;
	int j;

	if (NULL == rsa || NULL == pubkey || 0 == pubkey_size)
		return size;
	if (NULL != pubkey && pubkey_size >= size) {
		for (j = 0; j < RSALIMBS; j++)
			pkey[j] = htonl(r->n[j]);
		return 0;
	}
	errno = EINVAL;
	return -1;
}

/**
 * @brief Store a (third party) public key into a RSA context
 *
 * Store a (third party) public key into a RSA context. At
 * the same time reset the private key in the context, since
 * we don't have it.
 * If this function is called with either of rsa or pubkey NULL,
 * or with pubkey_size 0, the returned value is the number of
 * bytes wanted to complete a public key.
 * Otherwise the return value is zero, if the buffer had
 * the required number of bytes for the entire public key.
 * The return value is -1 and errno set to EINVAL if the buffer
 * was too small to keep a public key.
 *
 * @param rsa pointer to a RSA context allocated by rsa_init()
 * @param pubkey pointer to a buffer holding the public key
 * @param pubkey_size number of bytes in the buffer
 *
 * @result zero on success, size of public key when called with NULL, -1 error
 */
int rsa_set_public_key(void *rsa, void *pubkey, int pubkey_size)
{
	rsa_t *r = (rsa_t *)rsa;
	uint32_t *pkey = (uint32_t *)pubkey;
	int size = RSABITS / 8;
	int j;

	if (NULL == rsa || NULL == pubkey || 0 == pubkey_size)
		return size;
	if (NULL != pubkey && pubkey_size >= size) {
		for (j = 0; j < RSALIMBS; j++)
			r->n[j] = ntohl(pkey[j]);
		/* reset private key part */
		for (j = 0; j < RSALIMBS; j++)
			r->d[j] = 0;
		return 0;
	}
	errno = EINVAL;
	return -1;
}

/**
 * @brief Strip the topmost 2 bits (of RSABITS total) from a "plaintext"
 *
 * This is required in order to allow RSA encryption of a "plaintext",
 * because the modular exponentiation shall not overflow RSABITS.
 * The size must be a multiple of the RSALIMBS (= RSABITS/LIMBBITS),
 * and is in practise, in the case of Entropy, exactly RSALIMBS.
 * In theory you could encrypt multiple "plaintexts" at once by
 * applying the modular exponentiation to each block of RSALIMBS limbs.
 *
 * @param rsa pointer to a RSA context allocated by rsa_init()
 * @param plaintext pointer to a buffer to "strip"
 * @param size number of bytes in the buffer
 *
 * @result zero on succes, -1 and errno = EINVAL in case of wrong arguments
 */
int rsa_strip(void *rsa, void *plaintext, int size)
{
	rsa_t *r = (rsa_t *)rsa;
	limb_t *pt = (limb_t *)plaintext;
	int i;

	if (NULL == r) {
		errno = EINVAL;
		return -1;
	}
	/* reject odd sizes */
	if (0 != (size % (LIMBBITS/8))) {
		errno = EINVAL;
		return -1;
	}
	/* convert to number of limbs */
	size /= LIMBBITS / 8;

	/* reject non-multiples of RSALIMBS */
	if (0 != (size % RSALIMBS)) {
		errno = EINVAL;
		return -1;
	}
	for (i = 0; i < size; i += RSALIMBS) {
		limb_t msb = pt[i+RSALIMBS-1];
		pt[i+RSALIMBS-1] = msb & ~(DIGMSB | (DIGMSB >> 1));
	}
	return 0;
}

/**
 * @brief Decode a RSA ciphertext into a plaintext
 *
 * Decode the ciphertext into a plaintext, using the private key
 * from the given RSA context. The size must be a multiple of RSALIMBS.
 *
 * @param rsa pointer to a RSA context allocated by rsa_init()
 * @param ciphertext pointer to a RSA encrypted ciphertext
 * @param plaintext pointer to a buffer to receive the plaintext
 * @param size size of both, ciphertext and plaintext buffers, in bytes
 *
 * @result zero on success, -1 and errno = EINVAL in case of wrong arguments
 */
int rsa_decode(void *rsa, void *ciphertext, void *plaintext, int size)
{
	static limb_t pt[BN_MAXSIZE];
	static limb_t ct[BN_MAXSIZE];
	rsa_t *r = (rsa_t *)rsa;
	limb_t *psrc = (limb_t *)ciphertext;
	limb_t *pdst = (limb_t *)plaintext;
	int i, j;

	if (NULL == rsa) {
		errno = EINVAL;
		return -1;
	}

	if (0 != (size % RSALIMBS)) {
		errno = EINVAL;
		return -1;
	}

	if (NULL == ciphertext || NULL == plaintext) {
		return size;
	}

	for (i = 0; i < size; /* */) {
		for (j = 0; j < RSALIMBS; j++) {
			ct[j] = ntohl(*psrc++);
			i += sizeof(limb_t);
		}
		/* decrypt: pt = ct^d % n */
		bn_modexp(pt, ct, r->d, r->n, RSALIMBS);
		for (j = 0; j < RSALIMBS; j++)
			*pdst++ = pt[j];
	}
	memset(ct, 0, sizeof(ct));
	memset(pt, 0, sizeof(pt));
	r = NULL;
	psrc = NULL;
	pdst = NULL;
	return 0;
}

/**
 * @brief Encode a plaintext into a RSA ciphertext
 *
 * Encode the plaintext into a RSA ciphertext, using the public
 * key from the given RSA context. The size must be a multiple
 * of RSALIMBS.
 *
 * @param rsa pointer to a RSA context allocated by rsa_init()
 * @param plaintext pointer to a buffer with the plaintext
 * @param ciphertext pointer to buffer to receive the ciphertext
 * @param size size of both, plaintext and ciphertext buffers, in bytes
 *
 * @result zero on success, -1 and errno = EINVAL in case of wrong arguments
 */
int rsa_encode(void *rsa, void *plaintext, void *ciphertext, int size)
{
	static limb_t pt[BN_MAXSIZE];
	static limb_t ct[BN_MAXSIZE];
	rsa_t *r = (rsa_t *)rsa;
	limb_t *psrc = (limb_t *)plaintext;
	limb_t *pdst = (limb_t *)ciphertext;
	int i, j;

	if (NULL == rsa) {
		errno = EINVAL;
		return -1;
	}

	if (0 != (size % RSALIMBS)) {
		errno = EINVAL;
		return -1;
	}

	if (NULL == ciphertext || NULL == plaintext) {
		return size;
	}

	for (i = 0; i < size; /* */) {
		for (j = 0; j < RSALIMBS; j++) {
			pt[j] = *psrc++;
			i += sizeof(limb_t);
		}
		/* encrypt: ct = mi^e % n */
		bn_modexp(ct, pt, r->e, r->n, RSALIMBS);
		for (j = 0; j < RSALIMBS; j++)
			*pdst++ = htonl(ct[j]);
	}
	memset(pt, 0, sizeof(pt));
	memset(ct, 0, sizeof(ct));
	r = NULL;
	psrc = NULL;
	pdst = NULL;
	return 0;
}

#if	BN_TEST

int info(const char *fmt, ...)
{
	va_list ap;
	int len;

	va_start(ap, fmt),
	len = vfprintf(stdout, fmt, ap);
	va_end(ap);
	return len;
}

void tv_add(struct timeval *dst, const struct timeval *src, uint64_t add)
{
	uint32_t usec, sec;

	usec = (uint32_t)(add % 1000000);
	sec = (uint32_t)(add / 1000000);

	dst->tv_usec = src->tv_usec + usec;
	dst->tv_sec = src->tv_sec + sec;
	while (dst->tv_usec >= 1000000) {
		dst->tv_usec -= 1000000;
		dst->tv_sec += 1;
	}
}

void tv_sub(struct timeval *dst, const struct timeval *src, uint64_t sub)
{
	uint32_t usec, sec;

	usec = (uint32_t)(sub % 1000000);
	sec = (uint32_t)(sub / 1000000);

	dst->tv_usec = src->tv_usec - usec;
	dst->tv_sec = src->tv_sec - sec;
	while (dst->tv_usec < 0) {
		dst->tv_usec += 1000000;
		dst->tv_sec -= 1;
	}
}

void tv_div(struct timeval *dst, const struct timeval *src, uint32_t div)
{
	uint32_t rem;
	rem = src->tv_sec % div;
	dst->tv_sec = src->tv_sec / div;
	dst->tv_usec = src->tv_usec / div;
	dst->tv_usec += 1000000 * rem / div;
	while (dst->tv_usec > 1000000) {
		dst->tv_sec++;
		dst->tv_usec -= 1000000;
	}
}

#define	BN_ROUNDS	5

int bn_test(void)
{
	struct timeval tv0, tv1, tva;
	limb_t d[BN_MAXSIZE];
	limb_t e[BN_MAXSIZE];
	limb_t n[BN_MAXSIZE];
	limb_t mi[BN_MAXSIZE];	/* message in */
	limb_t mo[BN_MAXSIZE];	/* message out */
	limb_t ct[BN_MAXSIZE];	/* cipher text */
	cnt_t i, j;

	printf("Generating %d bits RSA keys\n", RSABITS);
	memset(&tva, 0, sizeof(tva));
	printf("RSA ");
	for (i = 1; i < BN_ROUNDS; i++) {
		printf("**** round #%d ****\n", i);
		gettimeofday(&tv0, NULL);

		bn_rsakey(n, e, d, 5);

		gettimeofday(&tv1, NULL);
		tv_sub(&tv1, &tv1, 
			(int64_t)1000000 * tv0.tv_sec + tv0.tv_usec);
		tv_add(&tva, &tva,
			(int64_t)1000000 * tv1.tv_sec + tv1.tv_usec);

		/* create a random message of RSABITS-2 (2 MSB are 0) */
		for (j = 0; j < RSALIMBS; j++)
			mi[j] = sl_random(0, LIMBMASK);
		mi[j-1] &= ~(DIGMSB | (DIGMSB >> 1));
		printf("mi: 0x%s\n", bn2x(mi, RSALIMBS));
		/* encrypt: ct = mi^e % n */

		bn_modexp(ct, mi, e, n, RSALIMBS);
		printf("ct: 0x%s\n", bn2x(ct, RSALIMBS));

		/* decrypt: mo = ct^d % n */
		bn_modexp(mo, ct, d, n, RSALIMBS);
		printf("mo: 0x%s\n", bn2x(mo, RSALIMBS));

		if (0 != bn_cmp(mi, mo, RSALIMBS)) {
			fprintf(stderr, "Decrypt error!\n");
			exit(1);
		}
	}
	printf("passed\n");
	tv_div(&tva, &tva, BN_ROUNDS);
	printf("avg. %ld.%06ld seconds\n", tva.tv_sec, tva.tv_usec);
	return 0;
}

int main(void)
{
	bn_test();
	return 0;
}
#endif
